An orbital- free molecular dynamics study of melting in K20, K55, K92, K142, Rbss and 

Cs 55 clusters. 
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The melting-like transition in potasium clusters Kjv, with N=20, 55, 92 and 142, is studied by 
using an orbital-free density-functional constant-energy molecular dynamics simulation method, 
and compared to previous theoretical results on the melting-like transition in sodium clusters of the 
same sizes. Melting in potasium and sodium clusters proceeds in a similar way: a surface melting 
stage develops upon heating before the homogeneous melting temperature is reached. Premelting 
effects are nevertheless more important and more easily established in potasium clusters, and the 
transition regions spread over temperature intervals which are wider than in the case of sodium. For 
all the sizes considered, the percentage melting temperature reduction when passing from Na to K 
clusters is substantially larger than in the bulk. Once those two materials have been compared for a 
number of different cluster sizes, we study the melting-like transition in Rbss and CS55 clusters and 
make a comparison with the melting behavior of Nass and K55. As the atomic number increases, 
the height of the specific heat peaks decreases, their width increases, and the melting temperature 
decreases as in bulk melting, but in a more pronounced way. 

PACS numbers: 36.40.Ei 64.70.Dv 



I. INTRODUCTION 

The melting transition of a pure bulk material occurs 
at a well defined temperature for a given external pres- 
sure. The melting-like transition in small clusters com- 
posed of a finite number of atoms spreads instead over a 
temperature interval that widens as the cluster size de- 
creases. A number of thermal properties sensitive to the 
cluster size and essentially different from those found in 
bulk materials emerge in the transition region defined by 
that temperature interval, which has motivated a lot of 
theoreticalErEl and experimentaJj~lL3 investigations. The 
experiments of Schmidt et a/Jl3 have shown strong non- 
monotonic variations of the melting temperature of free 
sodium clusters with size, which can not be completely 
explained either by electronic or geometric shell closing 
arguments. The theoretical simulations have predicted 
the occurrence of several premelting effects, like surface 
melting or structural isomerizations, and also the exis- 
tence of a dynamic coexistence regime, where the cluster 
can fluctuate in time between being completely solid or 
liquid. 

We have previously reported density functional 
orbital-free molecular dynamics (OFMD) simulations 
of the melting process Jju. sodium clusters Na_/v, pwith 
N=8,20, 55,92, and 142.LH3 The OFMD technique^ is 
completely analogous to the method devised by Car and 
Parrinello (CP-1 to perform dynamic simulations at an 
ab initio leveljlil but the electron density is taken as the 
dynamic-variable Jij as opposed to the Kohn-Sham (KS) 
orbitalsEEl in the original CP method. This technique, 
whose main advantage over KS-based methods is that 
the computational effort to update the electronic sys- 
tem scales linearly witfa-neluster size, IWfi-iJsefHi already 
used both in solid stateOlla and clusterB'Bflli3'E2l physics. 



Our predictions of the temperatures at which homoge- 
neous cluster melting occurs were in go od ag reement with 
the experiments of Haberland's group,t£lliil excepting the 
enhancement of the melting temperature around N=55, 
which was not reproduced. We also observed a number of 
interesting premelting effects, mostly the establishment 
of a surface melting stage at a temperature lower than the 
homogeneous melting temperature for Na2o, Nag2 and 
Nai42, and several isomerization transitions in Nas and 
Na2o- K is interesting to study similar systems like clus- 
ters of K, Rb and Cs in order to assess whether those 
trends are a general feature of alkali clusters or not. With 
this goal, we consider in this paper the melting-like tran- 
sition of Kn clusters, with N=20, 55, 92, and 142, and 
compare it with that of sodium clusters of the same size. 
As a second step, for a fixed cluster size of N=55, we 
study the melting behavior of Rb55 and CS55 and com- 
pare it with that of Na55 and K55 . In the next section we 
briefly present some technical details of the method. The 
results are presented and discussed in section III and, fi- 
nally, section IV summarizes our main conclusions. 



II. THEORY 

The orbital-free molecular dynamics method is a Car- 
Parrinello total energy schemeE-il which uses an explicit 
kinetic-energy functional of the electron density, and has 
the electron density as the dynamical variable, as op- 
posed to the KS single particle wavefunctions. The 
main features of the energy functional and the calcu- 
lation scheme-. haiie been described at length in previ- 
ous work,iZH3BE3 and details of our method are as 
described by Aguado et al.lD In brief, the electronic ki- 
netic energy functional of the electron density, n(f), cor- 
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responds to the gradient expansion, around the homoge- 
neous limit through second orderc3H~t3 



T s [n]=T lt [n] + -T w [n), 



(1) 



where the first term is the Thomas-Fermi functional 
(Hartrcc atomic units have been used) 



T TF [n} = ±(3nY /3 Jn(rf/*dr, (2) 

and the second is the lowest order gradient correction, 
where T w , the von Weizsacker term, is given by 



T w \n] = - 



Vn(r) 
n(r) 



-dr. 



(3) 



The local deasiky approximation is used for exchange and 
correlation.Bj'Ea In the external field acting on the elec- 
trons, V ext {r) = ^2 n v(r — R n ), we take v to be the lo- 
cal pseudopotential of Fiolhais et aL,c3 which reproduces 
well the properties of bulk alkalis and has been shown to 
have good transferability to alkali clusters.EH The clus- 
ter is placed in a unit cell of a cubic superlattice, and 
the set of plane waves periodic in the superlattice is used 
as a basis set to expand the valence density. Following 
Car and Parrinello,EJ the coefficients of that expansion 
are regarded as generalized coordinates of a set of ficti- 
tious classical particles, and the corresponding Lagrange 
equations of motion for the ions and the electron density 
distribution are solved as described in Ref. fj| 

The calculations used a supercell of edge 71 a.u. for 
K 2 o and 81 a.u. for K 55 , K 92 , Ki 42 , Rb 55 and Cs 55 . 
An energy cut-off of 8 Ryd was used in the plane wave 
expansion of the energy for K clusters, and of 6.15 Ryd 
for Rb and Cs clusters. In all cases, a 64x64x64 grid 
was used. The cut-offs used give a convergence of bond 
lengths and binding energies as good as that obtained 
for sodium clusters J3 The fictitious mass associated with 
the electron density coefficients ranged between 6.3xl0 8 
and 4.0xl0 9 a.u., depending on the material and on the 
temperature of the simulations. The equations of motion 
of K clusters were integrated using the Verlet algorithmo 
for both electrons and ions with a time step ranging from 
At = 0.83 x 10~ 15 sec. for the simulations performed at 
the lowest temperatures, to At = 0.67 x 10~ 15 sec. for 
those at the highest ones. In the case of Rb55 the time 
steps ranged from At = 2.38 x 10~ 15 sec. to At = 1.31 
x 10~ 15 sec, and in the case of CS55 from At = 4.29 
x 10~ 15 sec. to At = 1.79 x 10~ 15 sec. These choices 
resulted in a conservation of the total energy better than 
0.1 % in all cases. 

Several molecular dynamics simulation runs at differ- 
ent constant energies were performed in order to obtain 
the caloric curve for each cluster. The initial positions 
of the atoms for the first run were taken by slightly de- 
forming the equilibrium low-temperature geometry of the 
cluster. The final configuration of each run served as the 



starting geometry for the next run at a different energy. 
The initial velocities for every new run were obtained 
by scaling the final velocities of the preceding run. The 
total simulation times for each run at constant energy 
were 40 ps for K clusters. 80 ps for Rbs5, and 140 ps for 
CS55 . These different simulation times were chosen in or- 
der to obtain a good convergence of the several melting 
indicators described below. They increase with atomic 
number because the typical atomic vibrational frecuen- 
cies decrease with atomic number. 

A number of indicators to locate the melting-like tran- 
sition were employed. Namely, the specific heat defined 
byEl 



C v = [N — N(l - 



3N -6 



) < Ekin >t 



(4) 



where N is the number of atoms and <>* indicates the 
average along a traject«ry; the root-mean-square (rms) 
bond length fluctuationtl 



(< R% > t - < R^ >f) 1 /2 



N(N- 1) ^ 



< Rij >t 



the "atomic equivalence indexes" 



(5) 



(6) 



and finally, the average over a whole dynamical trajectory 
of the radial atomic distribution function g(r) , defined by 



dN at = g{r)dr 



(7) 



where dN a t(r) is the number of atoms at distances from 
the center of mass between r and r + dr. _ _ 

The experimental results of Haberland et a/.E2rE3 sug- 
gest that both electronic and atomic shell effects deter- 
mine the irregular size evolution of the melting temper- 
atures of sodium clusters. Our orbital- free method does 
not account for quantum shell effects, and gives cluster 
energies that vary smoothly as a function of cluster size 
(that is, it does not reproduce the energy oscillations as- 
sociated with electronic shell closures). All atomic shell 
effects associated with the geometrical arrangement of 
ions are properly accounted for, however. Thus, our 
method is not expected to give a detailed account of the 
size variation of the melting temperatures of metal clus- 
ters. It will work better for those cluster sizes where an 
electronic shell closing appears (N=20, 92 and 142), and 
worse for those sizes intermediate between two electronic 
shell closures (N=55). The material dependent trends 
for a given cluster size should be accurately predicted by 
our energy model, as these will be a consequence of the 
specific pseudopotential. 
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III. RESULTS 

A very important issue in the simulations of cluster 
melting is the election of the low-temperature isomer to 
be heated. A good knowledge of the ground state struc- 
ture (global minimum) is required, as the details of the. 
melting transition are known to be isomer-dependent.E2l 
But the problem of performing a realistic global opti- 
mization search is exponentially difficult as size increases, 
so finding the global minima of clusters with 55 atoms 
or more becomes impractical. In our previous works 
we directly started from icosahedral isomers forJ^ass, 
Na92 and pNa^, as there is some experimental] and 
theoretically indications that suggest icosahedral pack- 
ing in sodium clusters, and found a good agreement with 
the experimental results of Haberland's groupJid Simu- 
lated annealing runs for Nag2 and Nai42 always led to 
disordered structures with an energy higher than that of 
the corresponding icosahedral isomer. The melting be- 
havior of these disordered structures has been separately 
analyced ancLibund to be different from that of icosahe- 
dral clusters .lj As the comparison with experiment was 
favourable only for icosahedral isomers and the total en- 
ergy of these structures was always lower than that of 
disordered structures, we have chosen icosahedral isomers 
in the study of the melting behavior of large alkali clus- 
ters: K 55 , Rb 55 and Cs 55 are complete two-shell icosahe- 
drons, K92 and K442 are incomplete three-shell icosahe- 
drons constructed by following the icosahedral spewing 
pattern described by Montejano-Carrizales et a/£3 The 
low-temperature isomer of K20 was obtained by the dy- 
namic simulated annealing technique ,11-3 by heating the 
cluster to 400 K and then slowly reducing the tempera- 
ture. The resulting structure is essentially the same as 
that obtained for Na2o with the same technique!] 

The temperature evolutions of the specific heat C„ and 
of the rms bond length fluctuation 8 of K20 are shown in 
Fig 1. The specific heat displays two maxima around 90 
K and 130 K. The <5(T) curve has a small positive slope at 
low temperatures that reflects the thermal expansion of 
the solidlike cluster, and then two abrupt increases that 
correlate with the two peaks in the specific heat. Both 
magnitudes indicate that the melting of K20 occurs in two 
well separated steps over a wide range of temperatures. 
To analyce the nature of those two steps we show in Fig. 
2 short-time averages of the "atomic equivalence indexes" 
of K20 for a number of representative temperatures. For 
a temperature at which the cluster is completely solid, 
the <Tj(t) curves show a high degeneracy which is specific 
of the symmetry of the isomer under consideration!!!] The 
structure of K20, as that of Na2oJ3 can be divided into 
two subsets: two internal "core" atoms and 18 periph- 
eral "surface" atoms. The transition at k 90 K is iden- 
tified with an isomcrization transition in which the 18 
peripheral atoms begin to interchange their positions in 
the cluster while the two central atoms remain oscillat- 
ing around their initial positions. When a temperature 



of w 130 K is reached, one of the two inner atoms moves 
out to the cluster surface, while the other remains in 
its central position. Then the second transition is iden- 
tified with another isomcrization transition in which a 
new set of (19+1) isomers begins to be visited. The 8 
quantity increases with temperature after this point in a 
smooth way. This is due to the more and more frequent 
interchanges of the central atom with one of the periph- 
eral atoms upon increasing the temperature. Neverthe- 
less, the interchange rate between central and peripheral 
atoms remains slower than the interchange rate between 
peripheral atoms for all temperatures considered. 

Fig. 3 shows the specific heat and 8 curves for K55. 
The specific heat displays a main assymetric peak cen- 
tered approximately at 160 K, while 8 shows two abrupt 
increases at w 110 K and 160 K. The second abrupt in- 
crease in (5 coincides with the position of the main spe- 
cific heat peak. Although the first step in 8 is not in 
correspondence with any well-defined specific heat peak, 
there is a visible shoulder in the low temperature side (a 
clear assymetry) of that peak. Moreover, the width of 
the transition region, approximately 100 K, is predicted 
to be the same with both indicators. The nature of melt- 
ing is analyced by plotting the temperature evolution of 
the time-averaged radial atomic density distribution g(r). 
At a low temperature of T=86 K, Fig. 4 shows that the 
atoms are distributed in several icosahedral shells (in the 
outer shell, the twelve atoms in vertex positions can be 
distinguished from the rest due to slightly different radial 
distances). At the temperature where the first step in 5 
emerges, the detailed structure in the ionic density distri- 
bution has been washed out by the thermal effects, and 
the movies show that the cluster surface is melted. How- 
ever, the different shells are still clearly distinguished, 
showing that there are not interchanges between atoms 
in different shells. At a temperature higher than 160 K, 
the disctintion of the several radial shells is not possi- 
ble anymore. All the atoms are able to diffuse across 
the cluster volume, that is both intrashell and intershell 
displacements are allowed, and the liquid phase is com- 
pletely established. Upon a further increasing in temper- 
ature, the only appreciable change in g(r) is due to the 
thermal expansion of the cluster. 

The results for K92 are shown in Figs. 5 and 6. Both 
the specific heat and 8 predict a two-step melting pro- 
cess, with a first transition at ps 110 K and a second 
transition at w 200 K. As seen in Fig. 6, the first tran- 
sition is associated again with surface melting, with no 
substantial intershell diffusion. For temperatures higher 
than 200 K, the cluster is completely liquid. K142 melts 
in two main steps, at » 140 K and 230 K (Fig. 7). There 
is also a small bump on the low-temperature side of the 
first specific heat peak, correlated with a small abrupt 
increase of 8 at w 90 K. This previous step is associated 
with an isomerization regime in which different isomers 
preserving the icosahedral symmetry are visited, and was 
also found for Nai42.u These isomerizations involve the 
motion of the five vacancies in the outer shell. The sur- 
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face melting stage is not developed yet, however, as the 
icosahedral order persists. The distribution of atoms in 
three shells is still distinguished at a temperature of 190 
K where the cluster surface is melted. The average radial 
ionic density distribution is not completely uniform until 
the homogeneous melting temperature is reached- 

The bulk melting temperature of K (337 K)cJ is r£=. 
duced by a 10 % with respect to that of Na (371 K)H 
The melting temperaturea_pf K clusters are also smaller 
than those of Na clustersLra for all the sizes studied. The 
percentage reduction in melting temperature is substan- 
tially larger than in the bulk and a slightly decreasing 
function of cluster size (19 %, 17 %, 16 % and 15 % for 
N=20, 55, 92, and 142, respectively). The nature of the 
several premelting effects observed are similar for both 
materials. Nevertheless, they are more easily established 
in the case of K clusters. For example, surface and ho- 
mogeneous melting temperatures were closer together in 
the case of Nai42 (240 K and 270 K, respectively p than 
they are for K442 (140 K and 230 K, respectively); while 
two-step melting was not observed for a perfect two-shell 
Nass icosahedron,El the melting surface stage is well es- 
tablished for K55, which has the same low-temperature 
structure; in general, the transition region is wider for K 
than for Na clusters. 

As the main points in this comparison are quite inde- 
pendent of cluster size, we analyce in the following the 
melting behavior of Rb and Cs clusters for a fixed clus- 
ter size, namely N=55. Specific heat and S curves as a 
function of temperature are given in Figs. 9 and 10. The 
results are similar to those obtained for K55, namely a 
main assymetric specific heat peak and two steps in 5, the 
last of which correlates with the peak in the specific heat. 
The radial ionic density distributions of both materials 
present a similar temperature dependence and only the 
results for CS55 are shown in Fig. 11. The first transition 
(at w 110 K for Rb 55 and w 90 K for Cs 55 ) is identified 
with surface melting: Fig. 11 at 110 K shows that inter- 
shell diffusion is not important yet. The second (at « 140 
K for Rb55 and 130 K for CS55 ) corresponds to homoge- 
neous melting. We find that the different alkali clusters 
with N=55 atoms melt in a similar way. The main dif- 
ferences are the following: a) The homogeneous melting 
temperature decreases with increasing atomic number as 
in the bulk case, but in a more pronounced way. Specif- 
ically, in the series Na— >K^Rb^Cs, the bulk melting 
temperatures decrease byHsercentage values of 10 %, 7 
%, and 3 %, respectively,!^! while for the 55-atom clus- 
ters the corresponding percentage values are 17 %, 12.5 
%, and 7 %, respectively; b) The height of the specific 
heat peaks decreases and their width increases with in- 
creasing atomic number; c) Premelting effects are more 
important for the heavier elements. Specifically _two-step 
melting was not observed in the case of NassJa while a 
well defined surface melting stage is observed in the ther- 
mal evolution of K55, and CS55. 

It is perhaps not surprising that the melting tempera- 
ture reduction is larger in clusters compared to the bulk 



phase, where coordination effects associated with a large 
proportion of atoms in surface-like positions do not ap- 
pear. But a meaningful comparison can not be done due 
to the different structures adopted by alkali elements in 
the cluster (icosahedral packing) and bulk (bec packing) 
phases. The other two points do not invoke any com- 
parison with the bulk phase, r^pd can be more conve- 
niently addressed. Rey at alE3 have analyced the in- 
fluence of the softness of the repulsive core interaction 
on cluster melting. Specifically, a series of pair poten- 
tials differing just in their shape in the core region was 
constructed and used to investigate the melting behav- 
ior of 13-particle clusters. For those potentials with soft 
core repulsion, two-step melting was observed: the first 
step corresponds to the onset of frequent isomerizations 
involving only the surface atoms, while the second corre- 
sponds to homogeneous melting, involving also the cen- 
tral atom. For the harder potentials, those two steps 
merge into one, and melting-in-steps processes do not ap- 
pear. The repulsive part of our nscudopotential is harder 
the lighter is the alkali element ,0 so the importance of 
premelting effects can be expected to increase in the se- 
ries Na— >K— >Rb^Cs. In effect, a well-defined surface 
melting stage is not observed for Na55, while it develops 
before the homogeneous melting point-for the heavier al- 
kali elements. Moseler and NordiekEj have studied the 
influence of the potential range on the heat capacity of 
13-atom Morse clusters. They have found that decreas- 
ing the range of the potential increases the peak of the 
heat capacity in the melting transition region. The range 
of our pseudopotential increases with atomic number for 
the alkali elements.CJ We have performed a series of static 
calculations for the K2, Rb2 and CS2 molecules in order 
to construct their binding energy curves, and have found 
that the range of the interatomic interaction increases 
together with the range of the pseudopotentials, so a de- 
crease in the height of the specific heat peak is expected 
in the series Na— >K— >Rb— >Cs. This is what is observed 
indeed. Thus, we conclude that melting proceeds in a 
qualitatively similar way in clusters of the alkali elements 
Na, K, Rb and Cs, and that the small existing differences 
can be explained in terms of the different parameters 
defining the corresponding local pseudopotentials. 

A few comments regarding the quality of the simu- 
lations and of the annealing runs are perhaps in order 
here. The orbital-free representation of the atomic in- 
teractions is much more efficient than the more accurate 
KS treatments, but is still substantially more expensive 
computationally than a simulation using phenomenolog- 
ical many-body potentials. Such potentials contain a 
number of parameters that arc usually chosen by fitting 
some bulk and/or molecular properties. In contrast our 
model is free of external parameters, although there are 
approximations in the kinetic and exchange-correlation 
functionals. The orbital-free scheme accounts, albeit ap- 
proximately, for the effects of the detailed electronic dis- 
tribution on the total energy and the forces on the ions. 
We feel that this is particularly important in metallic 
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clusters for which a large proportion of atoms are on the 
surface and experience a very different electronic envi- 
ronment than an atom in the interior. Furthermore, the 
adjustment of the electronic structure and consequently 
the energy and forces to rearrangements of the ions is also 
taken into account. But the price to be paid for the more 
accurate description of the interactions is a less complete 
statistical sampling of the phase space. The simulation 
times are substantially shorter than those that can be 
achieved in phenomenological simulations. Nevertheless, 
we expect that the locations of the several transitions are 
reliable, because all the indicators we have used, both 
thermal and structural ones, are in essential agreement 
on the temperature at which the transitions start. 

IV. SUMMARY 

The melting-like transition in Kjy, with N=20, 55, 92 
and 142, Kb^ and CS55 clusters has been investigated 
by applying an orbital-free, density-functional molecu- 
lar dynamics method. The computational effort which 
is required is modest in comparison with the traditional 
Car-Parrinello Molecular Dynamics technique based on 
Kohn-Sham orbitals, that would be very costly for clus- 
ters of these sizes. The details of the several transitions 
have been explained and found to be similar to those 
found in the melting- like transition of sodium clusters flu 
alkali clusters show generally a separate surface melting 
stage prior to homogeneous melting. The homogeneous 
melting temperature has been found to decrease with in- 
creasing atomic number as in the bulk limit, but the per- 
centage value of this temperature reduction is larger than 
for the bulk materials, due to the larger proportion of 
surface-like atoms existing in the cluster; in fact, the re- 
duction in melting temperature when passing from Na to 
K clusters slightly decreases with increasing cluster size. 
The height of the specific heat peaks decreases and their 
width increases with increasing atomic number, and the 
premelting effects are more important for the heavier al- 
kalis. These trends have been rationalized in terms of 
physical features of the local pseudopotentials employed. 
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Captions of Figures. 

Figure 1 Specific heat (a) and S (b) curves of K20, tak- 
ing the internal cluster temperature as the independent 
variable. The deviation around the mean temperature is 
smaller than the size of the circles. 

Figure 2 Short-time averaged distances < 7"i(t) > s ta 



between each atom and the center of mass in K20, as 
functions of time for four representative values of the 
internal temperature. The bold lines follow the evolution 
of particular atoms. 

Figure 3 Specific heat (a) and 5 (b) curves of K55, tak- 
ing the internal cluster temperature as the independent 
variable. The deviation around the mean temperature is 
smaller than the size of the circles. 

Figure 4 Time averaged radial atomic density distri- 
bution of K55, at four representative temperatures. 

Figure 5 Specific heat (a) and 5 (b) curves of K92, tak- 
ing the internal cluster temperature as the independent 
variable. The deviation around the mean temperature is 
smaller than the size of the circles. 

Figure 6 Time averaged radial atomic density distri- 
bution of Kg2, at four representative temperatures. 

Figure 7 Specific heat (a) and 6 (b) curves of K142, 
taking the internal cluster temperature as the indepen- 
dent variable. The deviation around the mean tempera- 
ture is smaller than the size of the circles. 

Figure 8 Time averaged radial atomic density distri- 
bution of K142, at four representative temperatures. 

Figure 9 Specific heat (a) and 8 (b) curves of Rb^, 
taking the internal cluster temperature as the indepen- 
dent variable. The deviation around the mean tempera- 
ture is smaller than the size of the circles. 

Figure 10 Specific heat (a) and 8 (b) curves of CS55, 
taking the internal cluster temperature as the indepen- 
dent variable. The deviation around the mean tempera- 
ture is smaller than the size of the circles. 

Figure 11 Time averaged radial atomic density dis- 
tribution of CS55, at four representative temperatures. 
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